/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
/** @file Crout.h
 *  @brief Definition of class Crout - LU decomposition of general matrices
 *
 *  Crout's LU decomposition     A = L * U 
 *
 *  DIRECT SOLUTION OF LINEAR EQUATIONS A * x = b or A * X = B
 *  
 *  LU decomposition of general (non-symmetric) matrices used for:
 *  (1) Solution of linear equations
 *  (2) Inverse of a general matrix
 *  (3) Determinant of a general matrix
 * 
 */ 

#ifndef _CROUT_H
#define _CROUT_H

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <cassert>
#include <string>
#include "Arrays.h"
#include "iArrays.h"

using namespace std;

namespace smv_lib
{

//! small number - equal to zero
const double LU_TINY = 1.0e-20;

//---------------------------------------------------------------
/**
** \brief class Crout - LU decomposition of general matrices
**        Matrix A = L * U
**
**  DIRECT SOLUTION OF LINEAR EQUATIONS A * x = b
**
**  LU decomposition of general matrices is used for:
**  (1) Solution of linear equations
**  (2) Inverse of a general matrix
**  (3) Determinant of a general matrix
**  NOTE: Original matrix A is the SAME after LU decomposition
**/
//---------------------------------------------------------------
class Crout
{
  private:
    //! Matrix A (square real and positive definite)
    Matrix A;
    //! Order of matrix A
    int N;
    //! Vector b - right-hand side of equations A * x = b
    Vector b;
    //! Vector x - solution of linear equations A * x = b
    Vector x;
    //! Number of right-hand side vectors (no. of eqs. with the same coefficient matrix A)
    int M;
    //! Matrix B - right hand side of equations A * X = B (B is matrix N x M)
    Matrix B;
    //! Matrix X - the solution of linear equations A * X = B (X is matrix N x M)
    Matrix X;
    //! Inverse of the matrix A
    Matrix* A_inv;
    //! Determinant of the matrix A
    double Det;
    //! Lower and upper triangular factors (A = L * U)
    Matrix L, U;
    //! Indicator: true if matrix A is given
    bool isMatrix_A;
    //! Indicator: true if vector b is given (A*x=b)
    bool isVector_b;
    //! Indicator: true if matrix B is given (A*X=B)
    bool isMatrix_B;
    //! Indicator: true if more sistems of eqs. with the same matrix A are solved
    bool isMore;
    //! Indicator: true if decomposition is over
    bool isDecomposed;
    //! Indicator: true if inverse of the Matrix A is obtained
    bool isInverse;
    //! Indicator: true if linear equations are solved
    bool isSolved;
    //! output stream
    ofstream out;
    //! Name of the problem (output file is Name + ".txt")
    string Name;
    //! Indicator: true if the name of the problem is given
    bool isNameGiven;
    //! Working matrix - LU decomposition of the matrix A
    Matrix A_LU;
    //! Working vector in LU decomposition
    iVector LU_indx;
    //! Working double in LU decomposition
    double LU_d;

  public:
    //! empty (default) constructor
    Crout(void);
    //! constructor with a given matrix (to obtain LU decomposition only)
    Crout(const Matrix& AA);
    //! constructor with a given matrix A and vector b - to solve equation A*x=b
    Crout(const Matrix& AA, const Vector& bb);
    //! constructor with a given matrix A and matrix B - to solve equations A*X=B
    Crout(const Matrix& AA, const Matrix& BB);
    //! destructor
    ~Crout(void);

    //----------------------------- set and get functions
    //! set Matrix A (for empty constructor)
    void setMatrix_A(const Matrix& AA);
    //! set the name
    void setName(const string name) {
      Name = name;
      isNameGiven = true;
    }
    //! set the right-hand side vector b
    void setVector_b(const Vector& bb);
    //! set the right-hand side matrix B
    void setMatrix_B(const Matrix& BB);
    //! get lower triangular factor L (A = L * U)
    Matrix get_L();
    //! get upper triangular factor L (A = L * U)
    Matrix get_U();
    //! get inverse of the matrix A
    Matrix getInverse();
    //! get determinant of the matrix A
    const double getDeterminant();
    //! get solution of linear equations A * x = b
    Vector getSolution();
    //! get one solutions of the system of linear equations A * X = B
    Vector getSolutionNo(int i);
    //! get all solutions of the system of linear equations A * X = B
    Matrix getAllSolutions();

    //---------------------------------- working functions
    //! decomposition of the matrix A
    void Decomposition();
    //! inverse of the matrix A
    void Inverse();
    //! solution of linear equations A * x = b or A * X = B
    void Solution();
    //! write results to predefined file "CHOLESKY.txt"
    void Write();
    //! write results to file with a given name
    void WriteTo(const char* file);

  private:
    //! LU decompostition
    void ludcmp();
    //! LU backsubstitution
    void lubksb(Vector& bb);
    //! inverse matrix
    Matrix croutInv();
    //! solution of equations A * x = b
    Vector croutEqs(const Vector& bb);
    //! solution of a system of linear equations A * X = B
    Matrix croutEqs(const Matrix& BB);
    //! determinant of the matrix A
    const double Determinant();

};

//================ External functions for LU decopmposition
void LUdcmp(Matrix& A, iVector& indx, double& d);
void LUbksb(Matrix& A, iVector& indx, Vector& b);

}  // end of namespace smv_lib

#endif   // _CROUT_H

//================================ end of file
